Thermostatting by deterministic scattering 
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We present a mechanism for thermalizing a moving particle by microscopic deterministic scattering. 
As an example, we consider the periodic Lorentz gas. We modify the collision rules by including 
energy transfer between particle and scatterer such that the scatterer mimics a thermal reservoir 
with arbitrarily many degrees of freedom. The complete system is deterministic, time-reversible, 
and provides a microcanonical density in equilibrium. In the limit of the disk representing infinitely 
many degrees of freedom and by applying an electric field the system goes into a nonequilibrium 
steady state. 

PACS numbers: 05.20.-y, 05.45. +b, 05.60. +w, 05.70.Ln, 44.90.+C 



Thermostatting is a mechanism by which the internal 
energy of a many-particle system, and thus its temper- 
ature, are kept constant although there is a flux of en- 
ergy through the system as, e.g., induced by external 
fields, or by imposing temperature or velocity gradients 
0. In order to study nonequilibrium transport in fluids 
by computer simulations, Evans, Hoover, Nose and oth- 
ers developed methods of thermostatting by introducing 
a fictitious frictional force into the microscopic equations 
of motion [|]|| , modeling the interaction of particles with 
a thermal reservoir. This force is chosen to be velocity de- 
pendent in a way that the (total or kinetic) energy of the 
particles remains constant and that, in an equilibrium sit- 
uation, the system approaches microcanonical or canoni- 
cal distributions of the phase space variables. In contrast 
to stochastic thermostats (l],||||, this class of dynamical 
systems is completely deterministic and time-reversible, 
making them suitable models to study the connection 
between microscopic reversibility and macroscopic irre- 
versibility. This led to interesting new relations between 
statistical physics and dynamical systems theory P,|l^|, 
especially between transport coefficients and Lyapunov 
exponents j7j-|lC|], and between irreversible entropy pro- 
duction and phase space contraction P,^,|il]-|l3|] . On the 
other hand, employing these friction coefficients for creat- 
ing nonequilibrium steady states is physically somewhat 
obscure, because the equations of motion are fundamen- 
tally modified. It has been shown that the resulting dy- 
namical systems belong to some new class of generalized 
Hamiltonian systems but the problem still re- 

mains whether the relations mentioned above rely on hav- 
ing specifically this class of systems, or whether they are 
of a general nature ^,^-0. Much recent work 17 1^] 
has been devoted to deal with this question, and to relate 
this way to model nonequilibrium states to the one ini- 
tiated by Gaspard, where nonequilibrium is induced by 
appropriate boundary conditions in spatially extended 
systems @§0|. 

One of the simplest examples of a dynamical system 
in which thermostatting by a velocity-dependent fric- 
tion coefficient has been applied is the field driven pe- 
riodic Lorentz gas. The classical periodic Lorentz gas, 



where a point particle scatters elastically at hard disks 
arranged on a triangular lattice, has been the subject 
of many investigations and serves as a standard model 
in the field of chaos and transport |l6|-|20|]. In the field 
driven case an electric field acts on the moving particle 
and a thermostat keeps the energy of the particle con- 
stant at all times |^,|7 , |lO| , p^| , pTl -|24[| . For unit mass and 
unit charge of the particle the equations of motion read 
r = v , v = e — ((v)v , plus the geometric constraints 
induced by the disks, r and v represent position and ve- 
locity of the moving particle, e is the electric field, and 
stands for the friction coefficient which, by requiring 
energy conservation, is determined to be C( v ) = £ • v/v 2 . 
This set-up provides an example of a so-called Gaussian 
thermostat 

In this Letter we propose an alternative mechanism of 
deterministic, timc-revcrsiblc thermostatting that does 
not appeal to such a frictional force. We shall implement 
this mechanism by adapting the geometry of the periodic 
Lorentz gas, in which it suffices to consider one scatterer 
in an elementary cell supplemented by periodic bound- 
ary conditions, as depicted in Fig. |] (a). For the spacing 
between two neighboring disks at disk radius R = 1 we 
choose, following the literature @,|l^,||,|4) , w ~ 0.2361, 
ensuring that no particle can move collision-free for infi- 
nite time. Fig. [I] (b) defines the relevant variables for the 
collision process. At the collision we express the velocity 
v of the colliding particle in local polar coordinates 



(hi 





FIG. 1. (a) Elementary cell of the periodic Lorentz gas on 
a triangular lattice, (b) Definition of the relevant variables to 
describe the collision process. 
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(7,^), where 7 is the angle of incidence and v the abso- 
lute value of the velocity. The dashed variables indicate 
the respective values after the collision. We also intro- 
duce the angle (3, which determines the position of the 
colliding particle at the disk. For an elastic collision, as in 
case of the original Lorentz gas, one has (7, v) = (7', v'). 
In contrast to this, we propose to include an energy 
transfer between particle and disk, which, viewed in this 
way, serves as a momentum- and energy-reservoir. We 
do this by introducing an additional velocity variable 
k associated to the disk and by allowing that 7 ^ 7'. 
We require that the energy is conserved at the collision, 
v 2 + k 2 = v' 2 + k 12 = 2E, where E is the total energy 
of the system. Thus, the collision process in velocity 
space is still effectively defined by the dynamics of two 
variables for which we take (7, v). Using this setup, one 
can construct a model where the disk rotates with k as 
an angular velocity p5[ . By keeping the component of 
v perpendicular to the disk fixed and allowing only ex- 
change of energy via the tangent component, the col- 
lision process effectively reduces to the one of two col- 
liding masses on a line. Requiring energy and momen- 
tum conservation yields for the scattering rules a two- 
dimensional piecewise linear map. As a drastic simplifi- 
cation of such a model, but keeping important dynamical 
properties like time-reversibility, a deterministic dynam- 
ics, and the dynamical instability induced by the disk 
geometry, we choose here our collision rules according to 
a simple baker map |^J2^,^] . We apply it on the respec- 
tive Birkhoff coordinate of the ingoing angle, sin|7|, as 
its x-coordinate, and on v in the range of < v < \/2E 
as its y-coordinate. To ensure that the system is time- 
reversible, the forward baker acts if0<7<7r/2, and its 
inverse if — tt/2 < 7 < 0. The angle 7' always goes to 
the respective other side of the normal, as shown in Fig. 
|l| (b). For 7 > this gives 



M(sin|7|,u) 



(2sin| 7 |,t-/2) 



sin |t| < 0.5 



(2sin| 7 | - l,(v+ l)/2) ,sin| 7 | > 0.5 



(1) 



and vice versa for 7 < 0. As for k' , it is obtained from 
energy conservation. To avoid any symmetry breaking 
induced by this combination of forward and backward 
baker, we alternate their application in 7 with respect to 
the position f3 of the colliding particle on the circumfer- 
ence n. 

The above setting leads to a well-defined scattering sys- 
tem with three degrees of freedom. As it stands, however, 
it does not satisfy the microcanonical distribution, since 
the energy is not equipartitioned between all degrees of 
freedom, corresponding to a distribution p(v x ,v y , k) not 
being unifom on the energy shell v 2 + v 2 + k 2 = 2E. 
We incorporate this essential feature by amending the 
microscopic scattering rules, as given by the baker, in 
the most straightforward way: As a starting point, we 
calculate the projection of the microcanonical density 
p{v x ,v y ,k) — 5(2E — v 2 — v 2 — k 2 ) onto v yielding 



p(v) = vj y/2E(2E — v 2 ). We want that our system ap- 
proaches this density in the long time limit. Let p m ap(v) 
be the probability density for v at the moment of the col- 
lision corresponding to a respective Poincare surface of 
section, in contrast to the probability density of the time- 
continuous system, which we may denote with p con t(v)- 
During the free flight the particle cannot change its ve- 
locity, and thus p m ap{v) and p C ont(v) are simply related 
via the average time the particle travels between two col- 
lisions with velocity v. This average time plays the role 
of a weighting factor leading to 



Pcont(v) = -Pmapiy) 
V 



(2) 



where c is a constant to be fixed by normalization. Thus, 
we want that the map which determines the collision rules 
generates an invariant velocity distribution p* ont (v) = 
p(v) corresponding to 



P 



map 



(v) 



2v 2 



nEV2E - v 2 



(3) 



However, the invariant density of the baker map is sim- 
ply p*{x b ,Vb) = P*{xb)p*{vb) = 1- Therefore, we de- 
fine a conjugate map which produces the desired density 
by including a transformation Db = Y(v), where yg is 
the actual baker variable. This transformation must be 
continuous and invertible, as defined by conservation of 
probability [Eq], 



P*(VB)\dyB\ = Pma P (v)\dv\ . 
Y(v) can then immediately be computed to 

Y(V)- 



v J 2 v 

V 2E — v z H arcsm , 

ttE it J2E 



(4) 



(5) 



with < v < V2E , < Y(v) < 1. If we write x B = 
X(pf) = sin I7I, the full collision rules are thus given by 



( T ',»') = (r 1 ,r 1 )o¥o(i( 7 ),F(„)) 



(6) 



Computer simulations show that a Lorentz gas with these 
collision rules is microcanonical in both its position and 
momentum coordinates in phase space • 

We now inquire how the above ideas can be used to 
mimic the interaction of a moving particle with a thermal 
reservoir. For this purpose we associate arbitrarily many 
degrees of freedom to the disk which could be related, 
e.g., to different lattice modes in a crystal, as mechanisms 
for dissipating energy from a colliding particle. For sake 
of simplicity we do not distinguish here between all the 
individual velocities in the reservoir. Instead, we pretend 
that the particle interacts instantaneously with all (d—2) 
velocity components kj of the reservoir via an absolute 

reservoir velocity k = \Jk\ + k\ + . . . + k 2 _ 2 , to which 
we identify the disk velocity. We want that the projected 
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densities of the accessible variables (v x ,v y ,k) be gener- 
ated from the microcanonical distribution of the full d- 
dimensional system. In particular, the projection of the 
microcanonical distribution Pd(v x ,v y , k\, fe, • • • , fcd-2) = 
6(2E-vl-i 
for d > 2 to if 



J2j=i tf) onto pd(v) can be calculated 



Pd {v) = (d-2)(2E)- — v{2E- 



(7) 



Using the equipartition theorem E/d — T/2 with T the 
temperature and a Boltzmann constant ks — 1, and tak- 
ing the limit d — > 00, this expression reduces to the 
Maxwellian distribution Poo(v) — w/Texp(— v 2 l(2T)). 
Choosing p* ont . d (v) = pd{v) according to Eq. @ and 
using Eq. (|^) determines the corresponding density 
Pmap d( v ) °f the Poincare section. The transformation 
Yd(v) which yields p* map d (v) can then be calculated from 
Eq. (^) In the limit of d — > 00, Yoo^) becomes 

especially simple and reads 



(8) 



Y 00 (v) = -J-Z=ve -^ 2T '+erf 



with < v < 00 , < 3^(i>) < 1. Computer simula- 
tions show that a periodic Lorentz gas with these collision 
rules approaches projected densities in (v x ,v y ,k) which 
are identical to the ones obtained from a uniform dis- 
tribution on the energy shell in (v x , v y , kx, &2, • • • > 2) 
p5| . The temperature T of the equilibrium system is 
unambiguously defined via equipartitioning of energy as 
it enters into Eq. (||). Notice the complementarity be- 
tween our approach and the procedure adopted when us- 
ing stochastic boundary conditions P,p|,^0|. 

We are now ready to set up a nonequilibrium situa- 
tion by taking the system as defined in equilibrium for 
d — ► 00 and by switching on an electric field parallel to 
the x-axis. This field affects the velocity of the moving 
particle. However, since the particle is a small subsystem 
in a large reservoir, and since we have built in a mech- 
anism of equipartitioning of energy, the particle is still 
getting thermalized by our scattering rules with a tem- 
perature determined by the temperature T of the reser- 
voir. In fact, computer simulations show that the system 
approaches a nonequilibrium steady state with kinetic en- 
ergy and conductivity fluctuating around constant mean 
values ]25[ ] . That such a nonequilibrium steady state ex- 
ists according to our scattering mechanism is the central 
result of our Letter. 

In the following, we illustrate some important charac- 
teristics of this state by results obtained from computer 
simulations. Fig. ^ demonstrates that for small enough 
field strength and high enough temperature the energy of 
the system is still approximately equipartitioned between 
< v 2 > — < v x > 2 and the reservoir. Fig. || depicts the 
conductivity a(e) =< v x > je with respect to the field 
strength e. The strong decrease of tr(e) indicates that the 
system is in a highly nonlinear regime. Ohm's law may 
be suspected to hold only at very small field strengths 




FIG. 2. Relation between average velocity squared in the 



frame moving with the current, < v > 



< v x > 



and 



temperature T for the infinite dimensional model. Equiparti- 



tioning of energy would imply < v > 



< Vx > 



IT. 



e <C 0.1 |^2|,|TJ, however, in this regime reliable numeri- 
cal results are difficult to get. The broadest fluctuations 
on smaller scales are beyond numerical uncertainties and 
may be reminiscent to the strong irregularities as they 
occur in the conductivity of the Gaussian Lorentz gas 
|?],|l-24|. Unfortunately, it is not clear how to compare 
the conductivities of these two models quantitatively, 
since the choice of temperature in the Gaussian model 
is somewhat ambiguous by a factor of two (l2). Fig. |] 

(a) ,(b) show Poincare plots of (/3, sin.7) at the collisions. 
In (b) the deterministic baker has been replaced by a ran- 
dom number generator such that the system completely 
loses its memory at each collision. Fig. (a) indicates the 
existence of a fractal attractor, analogously to the one 
found in the Gaussian Lorentz gas |7|23 ; 3^], whereas in 

(b) the fractal structure is lost due to the stochastic- 
ity of the boundary conditions. Fig. ^ (c) should be 
compared to the analogous diagram obtained from the 
Lorentz gas with a Gaussian thermostat [p|, p2[p^ ]. In 
case of our model, there is no indication of a pruning- 
induced "bifurcation scenario" or an ergodic breakdown 
like in the Gaussian version. The same holds for other 
choices of Poincare sections in phase space. 

It would be interesting to compare our approach in 
more detail to the one based on frictional forces. For 
this purpose, a Lorentz gas with a so-called Nose-Hoover 
thermostat, where the particle velocity fluctuates around 
a mean value |Q, would be a more adequate model to 
compare with, but we are not aware that such a model 
has been studied. It would also be important to compute 
Lyapunov exponents and fractal dimensions of the phase 



FIG. 3. Conductivity cr(e) as it varies with field strength 
e. The curve consists of 90 data points, the numerical uncer- 
tainty of each point is less than 5 • 10 -4 . 
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(a) 



(b) 
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FIG. 4. (a),(b) Poincare section of (/3, sin 7) defined in Fig. 
|l]at the moment of the collision for field strength e = 1. In (a) 
the baker map of Eq. (jl[) has been used, in (b) the baker has 
been replaced by a random number generator, (c) Poincare 
section of /3 at the moment of the collision for varying field 
strength e. 

space densities for our model. This would enable to test 
the status of recently proposed connections between sta- 
tistical physics and dynamical systems quantities as dis- 
cussed in Refs. Jl(||J|l]|| . 

Our method has recently been applied to thermostat 
an interacting many-particle system of hard disks un- 
der shear j^], and with a temperature gradient, at the 
boundaries, where it leads to a connection between de- 
terministic thermostats and thermalization at stochastic 
boundaries [jjo). It would also be interesting to use it as 
a simple model for inelastic scattering in granular mate- 
rial, instead of employing velocity-dependent restitution 
coefficients Q. 

We are indebted to P.Gaspard and Chr. Wagner for 
many important hints. R.K. thanks the DFG for finan- 
cial support, K.R. thanks the European Commisssion for 
a TMR grant under contract no. ERBFMBICT96-1193. 
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